rsf_2022.R

Resource selection functions

Björn Reineking, 2022-09-09

Model the relative density of animals (also called range distribution or utilisation distribution) as a function of environmental predictors.

We will use the buffalo data set. Loading packages

library(animove)
library(ctmm)
library(sf)
library(mgcv)
library(mvtnorm)
library(lubridate)

Load buffalo data

See Kami’s slides for how we got here…

data(buffalo_utm)

Environmental data: topography, waterways, and NDVI

data(buffalo_env)
raster::plot(buffalo_env)

Animals and elevation

raster::plot(raster(buffalo_env, 1))
points(buffalo_utm)

To speed up analyses, we will only work with one individual, Cilla

cilla <- buffalo_utm[["Cilla"]]
raster::plot(raster(buffalo_env, 1), ext = extent(cilla) * 2)
lines(cilla)

The first day of Cilla shows some unrealistic movements. We remove those observations

cilla <- cilla[timestamps(cilla) > min(timestamps(cilla)) + days(1), ]

raster::plot(raster(buffalo_env, 1), ext = extent(cilla) * 2)
lines(cilla)

Minimal example of rsf.fit

Create telemetry object

cilla_telemetry <- as.telemetry(cilla)

## Warning in spTransform(xSP, CRSobj, ...): NULL source CRS comment, falling back
## to PROJ string

## Warning in wkt(obj): CRS object has no comment

## Warning in wkt(obj): CRS object has no comment

## Minimum sampling interval of 3 minutes in Cilla

Fit ctmm model

cilla_guess <- ctmm.guess(cilla_telemetry, CTMM=ctmm(isotropic = TRUE), interactive = FALSE)
cilla_select <- ctmm.select(cilla_telemetry, cilla_guess)

Fit akde

cilla_akde <- akde(cilla_telemetry, cilla_select)
plot(cilla_akde)

Create named list of rasters We can do this by hand

be <- list("elev" = raster(buffalo_env, "elev"),
            "slope" = raster(buffalo_env, "slope"),
            "var_NDVI" = raster(buffalo_env, "var_NDVI"))

The integrator = “Riemann” option is still experimental, we use it here because it is much faster

cilla_rsf_riemann <- rsf.fit(cilla_telemetry, cilla_akde, R = be, integrator = "Riemann")

## Warning in .local(x, ...): This function is only useful for Raster* objects with
## a longitude/latitude coordinates

## Maximizing likelihood @ n=166037

## Calculating Hessian

summary(cilla_rsf_riemann)

## $name
## [1] "OUF"
## 
## $DOF
##      mean      area diffusion     speed 
##  8.322854 16.577732 17.688692 17.345680 
## 
## $CI
##                                             low          est         high
## var_NDVI (1/var_NDVI)             -188.04892994 -72.68984716  42.66923562
## slope (1/slope)                     -0.37986892  -0.03499126   0.30988641
## elev (1/elev)                       -0.04693589  -0.01277037   0.02139514
## area (square kilometers)           230.24198135 398.32893545 611.73047276
## τ[position] (days)                   4.56389303   7.82598060  13.41967744
## τ[velocity] (minutes)               40.67981533  43.18035421  45.83459818
## speed (kilometers/day)              10.35589660  13.52672089  16.69123828
## diffusion (square kilometers/day)    3.12109675   5.29364725   8.03075027

A suitability map

suitability_riemann <- suitability(cilla_rsf_riemann, be, crop(be[[1]], extent(cilla) * 2))
raster::plot(suitability_riemann)

Range distribution (includes the ranging behaviour)

agde_cilla <- agde(cilla_rsf_riemann, be)
plot(agde_cilla)

Predict function

rsf.predict <- function(model, object, include_avail = TRUE) {
  log_avail <- 0
  if(include_avail) {
    xy <- as.data.frame(xyFromCell(object, 1:ncell(object)))
    xy <- sf::st_as_sf(xy, coords = c("x", "y"), crs = crs(object))
    xy <- st_transform(xy, CRS(model@info$projection))
    xy <- st_coordinates(xy)
    xy <- cbind(xy, -(xy[,1]^2 + xy[,2]^2)/2)
    beta_rr <- 1/model$sigma[1,1]
    beta_xyr <- c(model$mu*beta_rr, beta_rr)
    log_avail <- xy %*% beta_xyr
  }
  newdata <- as.data.frame(getValues(object))
  attr(newdata, "na.action") <- "na.pass"
  X <- stats::model.matrix.default(model$formula, data = newdata, na.action = na.pass)
  lambda <- X[,-1] %*% model$beta + log_avail # model.matrix.default puts as first column (Intercept) if model formula does not include "-1"
  r <- raster(object, 1)
  r[] <- exp(lambda - max(lambda, na.rm = TRUE))
  r <- r / sum(getValues(r), na.rm = TRUE)
  r
}

cilla_predict <- rsf.predict(cilla_rsf_riemann, crop(buffalo_env, extent(cilla) * 2), include_avail = FALSE)
raster::plot(cilla_predict)

Traditional RSF with downweighted Poisson regression

Functions to generate quadrature points and predict with the model

rsf_points <- function(x, UD, R = NULL, n = 1e5, k = 1e6, type = "Riemann", 
                       rmax =6*sqrt(UD@CTMM$sigma[1,1]),
                       interpolation = FALSE) {
  # Samples background points from a 2D normal distribution fitted to the relocation data, and extracts environmental
  # information from a raster object "R"
  # x: telemetry object
  # UD: UD object
  # R: raster* object
  # n: number of background points to sample
  # k: weight of presence points
  # rmax: maximum distance for Riemann-type integration
  # interpolation: do interpolation when sampling the grid
  # When type 0´= "MonteCarlo", importance sampling is done
  stopifnot(UD@CTMM$isotropic)
  stopifnot(type %in% c("Riemann", "MonteCarlo"))
  if (type == "Riemann") {
    quadrature_pts <- getValues(R)
    xy <- as.data.frame(xyFromCell(R, 1:ncell(R)))
    xy <- sf::st_as_sf(xy, coords = c("x", "y"), crs = crs(R))
    xy <- st_transform(xy, CRS(UD@info$projection))
    xy <- st_coordinates(xy)
    r <- sqrt(((xy[,1] - UD@CTMM$mu[1]))^2 + ((xy[,2] - UD@CTMM$mu[2]))^2)
    bg <- data.frame(case_ = 0,
                     x_ = xy[r<rmax,1], y_ = xy[r<rmax,2],  w_ = prod(res(R)), k_ = k)
    bg <- cbind(bg, quadrature_pts[r<rmax,])
    bg <- sf::st_as_sf(bg, coords = c("x_", "y_"), crs = CRS(UD@info$projection))
    xx <- data.frame(case_ = 1, x_ = x$x, y_ = x$y,
                     w_ = 1/k * UD$weights * mean(UD$DOF.area), k_ = k)
    xx <- sf::st_as_sf(xx, coords = c("x_", "y_"), crs = CRS(UD@info$projection))
    xx[names(R)] <- as.data.frame(raster::extract(R, st_transform(xx, crs(R)), method = ifelse(interpolation, "bilinear", "simple")))
    xx <- rbind(bg, xx)
  } else {
    quadrature_pts <- MASS::mvrnorm(n, mu = UD@CTMM$mu, Sigma = UD@CTMM$sigma)
    xx <- data.frame(case_ = 0, x_ = quadrature_pts[, 1], y_ = quadrature_pts[, 2], w_ = UD@CTMM$sigma[1,1]/n, k_ = k)
    xx <- rbind(xx, data.frame(case_ = 1, x_ = x$x, y_ = x$y,
                               w_ = 1/k * UD$weights * mean(UD$DOF.area), k_ = k 
    ))
    xx <- sf::st_as_sf(xx, coords = c("x_", "y_"), crs = CRS(UD@info$projection))
    xx[names(R)] <- as.data.frame(raster::extract(R, st_transform(xx, crs(R)), method = ifelse(interpolation, "bilinear", "simple")))
  }
  xy <- st_coordinates(xx)
  colnames(xy) <- c("x_", "y_")
  sd <- sqrt(UD@CTMM$sigma[1,1])
  xy[,1] <- (xy[,1] - UD@CTMM$mu[1])/sd
  xy[,2] <- (xy[,2] - UD@CTMM$mu[2])/sd
  xx <- cbind(xy, xx)
  xx
}

predict_rsf <- function(model, UD, object, include_avail = TRUE, data_crs = UD@info$projection) {
  if(include_avail) {
    xy <- as.data.frame(xyFromCell(object, 1:ncell(object)))
    xy <- sf::st_as_sf(xy, coords = c("x", "y"), crs = crs(object))
    xy <- st_transform(xy, data_crs)
    xy <- st_coordinates(xy)
    colnames(xy) <- c("x_", "y_")
    
    sd <- sqrt(UD@CTMM$sigma[1,1])
    xy[,1] <- (xy[,1] - UD@CTMM$mu[1])/sd
    xy[,2] <- (xy[,2] - UD@CTMM$mu[2])/sd
    
  } else {
    xy <- cbind("x_" = rep(0, ncell(object)), "y_" = rep(0, ncell(object)))
  }
  newdata <- as.data.frame(cbind(xy, getValues(object)))
  lambda <- as.numeric(predict(model, newdata, type = "link"))
  r <- raster(object, 1)
  r[] <- exp(lambda - max(lambda, na.rm = TRUE))
  r <- r / sum(getValues(r), na.rm = TRUE)
  r
}

A minimal “classic” example

Generate quadrature points (“background points”)

set.seed(2)
rsf_cilla_df <- rsf_points(cilla_telemetry, cilla_akde, buffalo_env, interpolation = TRUE)
rsf_cilla_df <- rsf_cilla_df[!is.na(rsf_cilla_df$slope),]

Fit a downweighted Poisson regression

m_rsf_cilla <- glm(case_*k_ ~ x_ + y_ + I(-(x_^2 + y_^2)/2) + elev + slope + var_NDVI, 
                  family = poisson(), data= rsf_cilla_df, weights = w_)

## Warning: glm.fit: fitted rates numerically 0 occurred

Summary of model and confidence intervals for parameter estimates

summary(m_rsf_cilla)

## 
## Call:
## glm(formula = case_ * k_ ~ x_ + y_ + I(-(x_^2 + y_^2)/2) + elev + 
##     slope + var_NDVI, family = poisson(), data = rsf_cilla_df, 
##     weights = w_)
## 
## Deviance Residuals: 
##      Min        1Q    Median        3Q       Max  
## -0.06440  -0.00419  -0.00044  -0.00005   0.59143  
## 
## Coefficients:
##                      Estimate Std. Error z value Pr(>|z|)    
## (Intercept)         -10.32153    3.56120  -2.898  0.00375 ** 
## x_                    0.31218    0.36103   0.865  0.38722    
## y_                    0.48704    0.44872   1.085  0.27775    
## I(-(x_^2 + y_^2)/2)   1.14409    0.27988   4.088 4.35e-05 ***
## elev                 -0.01275    0.01742  -0.732  0.46426    
## slope                -0.03507    0.17592  -0.199  0.84201    
## var_NDVI            -72.71439   58.86148  -1.235  0.21670    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 1125.8  on 289061  degrees of freedom
## Residual deviance: 1057.4  on 289055  degrees of freedom
## AIC: 1071.4
## 
## Number of Fisher Scoring iterations: 24

confint(m_rsf_cilla)

##                             2.5 %      97.5 %
## (Intercept)          -16.63420393 -2.92139831
## x_                    -0.38125884  1.06078566
## y_                    -0.33463784  1.41786282
## I(-(x_^2 + y_^2)/2)    0.67729735  1.77733492
## elev                  -0.04876778  0.01689897
## slope                 -0.48964633  0.20232405
## var_NDVI            -185.72957712 45.85066051

Map of suitability, including the home ranging behaviour

suitability_glm <- predict_rsf(m_rsf_cilla, cilla_akde, crop(buffalo_env, extent(cilla) * 2))
raster::plot(suitability_glm)

Map of suitability, without the home ranging behaviour

suitability_no_avail_glm <- predict_rsf(m_rsf_cilla, cilla_akde, crop(buffalo_env, extent(cilla) * 2), include_avail = FALSE)
raster::plot(suitability_no_avail_glm)

Multiple animals

• with rsf.fit: you can use the mean function on a list of rsf.fit objects
• “classic” approach: use glmmTMB and a mixed-effects Poisson regression: Muff, Signer & Fieberg (2020) J Anim Ecol 89: 80-92.